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0.  Introduction  and  Accomplishments 
Introduction 

■ -^The  recent  appearance  of  new  observational  data  from  specialized  satellites 
and  rocket  probes  has  led  to  increased  interest  in  upper  atmospheric  processes. 
The  work  to  be  reported  herein  is  of  the  current  status  of  a  limited  three- 
dimensional  model  of  the  dynamical  and  important  chemical  processes  which  are 
known  to  take  place  in  the  mesosphere  and  lower  thermosphere  to  an  altitude  of 
about  400  km  above  the  earth's  surfaces-  Unfortunately ,  funds  for  the  program 
were  cut  off  prior  to  its  completion  and  thus  the  model  codes  have  not  been 
finalized.  It  is  hoped  that  this  program  can  be  picked  up  again  in  the  near 
future. 

*^>The  modeling  approach  taken  was  to  make  use  of  the  dynamical  schemes  and 
simplified  chemical  treatments  embodied  in  our  three-dimensional  Stratospheric 
Circulation  Model  (SCM)  developed  for  the  study  of  stratospheric  ozone#'(see 
^unnold,  et  al_. ,  1975).  This  model  has  been  running  on  the  now  defunct  ILLIAC-4 


vector  computer  at  NASA's  Ames  Research  Center  in  California.  In  addition  to 
large  changes  required  in  the  existing  dynamics  and  chemistry  to  reform  the 
model  for  thermospheric  and  mesospheric  levels,  it  was  also  necessary  to  revise 
the  code  structure  to  accommodate  the  shift  from  the  ILLIAC  machine  to  the  AFGL 
CDC-660  computer.  Much  of  this  work  was  accomplished  on  the  NASA  machines  prior 
to  the  availability,  through  special  modems  and  telephone  connections,  of  the 
AFGL  CDC-6600.  Since  that  time  the  programs  have  been  transferred  to  AFGL  and, 
while  not  completed,  tests  of  the  model  dynamics  on  that  machine  have  been  under¬ 
taken.  —  _ _ _ 

The  basic  strategy  in  the  modeling  effort  was  to  use  the  modified  SCM  codes 
to  specify  the  large  scale  dynamical  properties  of  the  upper  atmospheric  region 
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and  for  the  integration  of  the  time  dependent,  three-dimensional  mass  continuity 
equations  for  the  chemically  active  species.  Development  of  the  chemical  and 
sub-scale  transport  properties  of  the  model  were  to  be  undertaken  by  the  AFGL 
group  under  the  direction  of  Dr.  S.  P.  Zimmerman.  Thus,  the  complete  modeling 
program  was  devised  to  be  a  cooperative  venture  between  Georgia  Tech  and  AFGL. 
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Accomplishments. 

The  program's  goal  was  to  create  a  three-dimensional  model  of  the  meso¬ 
sphere  and  lower  thermosphere  over  a  three-year  period  with  limited  funds.  The 
model  was  to  incorporate  simplified  dynamics  and  interactive  chemistry.  A 
"first"  run  of  a  single  simulation  experiment  with  the  completed  model  was  antic¬ 
ipated  late  in  the  third  year.  Thus,  intermediate  results  of  a  scientifically 
viable  nature  could  not  be  expected  prior  to  completion  of  the  program.  While 
the  program  has  been  cut  short  before  these  goals  could  be  attained,  substantial 
progress  has  been  made,  particularly  in  the  modification  of  the  dynamical  por¬ 
tions  of  the  model  codes  and  the  changes  required  to  run  the  model  on  the  AFGL 
CDC-6600  machine. 

A.  Model  heating 

One  of  the .principal  needs  of  the  upper  atmospheric  model  is  the  incorpora¬ 
tion  of  realistic  heat  forcing  processes  for  the  40-400  km  regions  of  interest. 
Considerable  thought,  therefore,  has  been  given  to  this  problem. 

The  existing  Dynamical/Chemical  Stratospheric  Circulation  Model  (SCM), 
which  is  being  revised  for  the  present  work  to  include  mesospheric  and  lower 
thermospheric  levels,  currently  uses  heating  codes  applicable  to  altitudes  be¬ 
low  %  80  km.  Since  a  number  of  quite  different  physical  processes  lead  to  atmo¬ 
spheric  heating  in  the  thermosphere,  new  model  codes  will  have  to  be  developed 
for  these  thermospheric  levels.  For  this,  however,  we  must  keep  in  mind  that 
such  codes  must  be  considerably  simplified  because  of  time  and  size  limitations 
which  must  be  imposed  on  the  already  large  three-dimensional  calculations. 

Some  of  the  heating  processes  we  have  considered  include: 
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(1 )  Direct  solar  absorption 

It  is  convenient  to  divide  the  solar  spectrum  into  several  large  wavelength 
segments  according  to  their  absorption  characteristics  and  treat  each  of  these 
segments  separately. 

a.  2050  -  3000  A  region. 

Heating  is  due  to  absorption  by  both  O2  and  0^  in  this  region  and 
is  most  important  at  the  lowest  thermospheric  and  mesospheric  levels 
(i.e.,  <  100  km).  Model  treatment  for  this  region  can  be  essentially 
the  same  as  was  devised  for  the  SCM  (see  Cunnold,  et^  aK  ,  1975).  For 
example,  for  03,  the  rate  of  temperature  change  due  to  0^  absorption  is 
approximated  by  the  linear  law 

X0 

^3t^03=  ~M~  Q(Nsec^) 

where  Xq3  is  the  number  mixing  ratio  of  03>  M  the  mass  of  an  average  air 
molecule,  Q(Nsec^)  the  heating  rate  due  to  absorption  by  one  molecule  of 
03,  N  is  the  number  of  03  molecules  in  the  cm  vertical  column  above  the 
point  of  heating,  and  c,  is  the  solar  zenith  angle.  In  the  SCM  calcula¬ 
tions  the  heating  rate,  Q,  is  approximated  by  a  finite  sum  over  a  number 
of  small  spectral  intervals  centered  on  wavelengths  in  the  form 

Q(Nsecc)  =  2  Oq3(A^ )  F(A,. )  exp(-  Nsecc) 

in  which  cxq^ ( )  is  the  absorption  coefficient  of  03  and  F( )  is  the 
solar  flux  of  photons  integrated  over  each  A^  interval.  In  the  calcu¬ 
lations,  tables  of  cig^,  F,  and  the  exponential  functions  are  maintained 
for  a  wide  range  of  likely  values. 
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b.  1027  -  1300  A  and  1750  -  2050  A  (Schumann-Runge  bands)  regions. 

Absorption  by  O2  in  these  two  banded  regions  is  comparatively  large 
below  about  120  km  altitude  but  can  probably  be  safely  neglected  above 

O 

this  level.  As  in  the  2050  -  3000  A  region  we  can  estimate  the  heating 
rates  by  summing  over  the  important  absorption  bands  using  band  averages 
as  tabulated  by  Hudson  and  Mahle  (1972)  for  the  1750  -  2050  A  region 
(although  the  apparent  temperature  dependence  of  the  cross  sections  is  a 

o 

complication)  and  by  Adams  (1974)  for  the  1027  -  1300  A  region. 

O 

c.  1300  -  1750  A  (Schumann-Runge  continuum)  regions. 

For  this  0^  absorption  region,  the  absorption  cross  sections  are 
quite  consistent  and  we  should  be  able  to  treat  this  region  as  a  single 
band.  Heating  by  absorption  in  this  region  is  most  important  to  total 
heating  at  levels  between  n,  100  and  130  km. 

d.  40  -  1027  A  ( EUV )  region. 

Nearly  all  photons  in  this  frequency  range  are  absorbed  by  photo¬ 
ionization  of  ,  O2,  and  0  which  leads  to  very  complicated  ionization 
and  photoelectric  processes.  These  processes  are  particularly  dominant 
above  n,  110  km  but  are  replaced  by  heating  through  collisional  processes 
above  0,  300  km.  It  may  be  possible  to  estimate  the  magnitude  of  the 

o 

heating  resulting  from  photon  absorption  in  the  40  -  1027  A  region  by 
using  a  simple  electron  density  model  such  as  that  of  Ching  and  Chiu  (1973) 
to  infer  photon  absorption  quantities  and  apply  a  heating  efficiency  fac¬ 
tor  of  n,  30  -  3511  (Stolarski,  e_t  al_. ,  1  975).  Such  a  model  is  currently 
being  tested. 
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(2 )  Atomic  oxygen  recombination  and  deactivation 

Atomic  oxygen  produced  at  high  model  altitudes  does  not  recombine  (and 
thus  release  its  chemical  energy  of  recombination)  above  n,  120  km.  Two  pro¬ 
cesses,  0  +  0+  M-*-02+M  and  0  +  02+M-*-02+M  may  be  important  here.  A 
simple  estimate  for  heating  by  these  processes  may  be  possible  by  assuming 
(Adams,  1974)  that  the  lifetime  of  an  oxygen  atom  goes  from  n,  5000  years  at 
150  km  to  %  2  hours  at  80  km.  Clearly,  model  vertical  transports  will  play  a 
large  role  here. 

The  process  of  deactivation  of  0(^D)  is  somewhat  uncertain  but  potent  . ;y 
important  to  heating  in  the  lower  thermosphere.  Whatever  the  mechanism 
(0(^0)  +  M  -*■  0(3P)  +  M  +  KE  is  the  prime  candidate),  the  reaction  takes  place 
very  fast  and,  since  there  is  no  known  large  source  for  0(^D)  at  night,  is 
confined  to  sunlit  hours.  A  possible  estimate  for  0(^D)  deactivation  in  the 
model  may  be  obtained  by  using  the  results  of  Adams  (1974,  pg.  97)  with  suit¬ 
able  adjustments  for  diurnal  and  latitudinal  variations. 

3 .  Molecular  thermal  conduction 

The  flux  of  heat  across  a  horizontal  surface,  Fz,  is  usually  parameter¬ 
ized  using 


where  here  K  represents  a  thermal  conduction  coefficient.  This  is  a  fairly 
simple  process  to  represent  computationally  but  the  selection  of  the  proper 
K's  will  be  done  in  consultation  with  the  AFGL  group. 

4.  15u  COq  and  62u  0  radiational  cooling 

Various  authors  have  estimated  cooling  rates  for  these  two  frequencies 
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and  the  model  parameterization  will  make  use  of  a  simplified  version  of  one  of 
these.  The  62 p  0  band  is  particularly  effective  above  %  110  km  while  the  15p 
CO^  band  seems  to  dominate  below  that  level. 

5.  Other  heating  processes 

Thermal  heating  by  the  dissipation  of  tidal  and  gravity  waves  may  be  impor¬ 
tant  to  the  thermosphere.  In  three  dimensional  models,  such  as  we  are  working 
with  here,  these  heat  quantities  will  be  realized  through  the  thermal  and  dynamic 
dissipation  terms  built  into  the  model  equations. 

Joule  heating  and  non-thermal  emissions  may  also  play  a  role  in  total  therm¬ 
ospheric  heating  but  their  magnitudes  are  usually  thought  to  be  small  and  will 
thus  be  neglected  in  the  current  work. 


B.  Model  lower  boundary  conditions 


The  Stratospheric  Circulation  Model  (SCM)  has  been  reconfigured  in  its  ver¬ 
tical  structure  to  incorporate  the  region  n.  40  km  -  400  km  in  its  26  vertical 

levels.  As  required  by  the  model,  mean  global  temperatures  (T)  and  stability 
dT  R  - 

quantities  -  q-  T)  at  each  level  were  obtained  from  the  U.S.  Standard  Atmo¬ 
sphere.  The  quantities  will  be  discussed  and  displayed  in  some  detail  in  later 
sections.  However,  we  want  to  point  out  that  the  new  Mesospheric  and  Lower 
Thermospheric  Model  (MTM)  overlaps  the  height  range  of  the  SCM  over  the  MTM's 
lowest  six  levels.  Thus,  it  will  be  possible  to  "drive"  the  lower  boundary  of 
the  MTM  using  values  computed  from  annual  runs  of  the  SCM.  To  this  end,  a  spec¬ 
ial  run  of  the  SCM  for  a  two  year  integration  period  was  performed  on  the  ma¬ 
chine  at  NASA's  Ames  Research  Center.  From  these  results,  we  have  obtained  for 
transference  to  AFGL: 

(1)  A  set  of  lower  boundary  conditions  for  temperature  (T),  vertical 
motion  (W)  and  ozone  (x).  A  complete  one-year  cycle  of  these  quantities  for 
the  model's  70  horizontal  degrees  of  freedom  were  collected  at  four-hour  inter¬ 
vals.  Thus,  we  have  tabulated  (on  a  computer  tape)  the  required  lower  boundary 
conditions  to  drive  the  MTM  as  functions  of  both  time  and  space.  This  involves 
more  than  1/2  million  values. 

(2)  Twelve  sets  of  initial  conditions,  one  for  each  month  of  the  year, 
were  generated  by  the  SCM  runs  and  tabulated  on  tape  files.  This  data  includes 
values  for  the  model  temperatures ,  vertical  motions  and  ozone  mixing  ratios  in 
the  region  of  overlap  between  the  SCM  and  the  MTM.  In  addition,  a  set  of  time 
dependent  total  heating  values  from  the  SCM  have  been  collected  for  the  one  year 
cycle  for  use  in  driving  the  MTM  during  early  dynamical  tests.  These  functions 
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will  be  replaced  by  internally  derived  heating  quantities  in  the  MTM's  final 
form. 

All  of  the  data  fields  described  in  (1)  ana  ^2)  have  been  transferred  (in 
ASCII  codes)  to  the  AFGL  6600  disk  system  and  are  available  for  use  in  the  model 
although  some  may  not  have  been  rewritten  in  binary  form  as  required  by  the  MTM 
input  scheme. 

To  incorporate  these  lower  boundary  conditions,  the  MTM  codes  have  been 
extensively  rewritten  and  tested.  Furthermore,  new  codes  have  been  generated  to 
allow  for  the  introduction  of  additional  minor  species  into  the  model  calcula¬ 
tion  in  fully  predictive  form  (through  the  species  continuity  equations).  De¬ 
tails  of  the  chemical  production  and  loss  terms,  however,  are  to  be  added  later 
in  cooperation  with  the  AFGL  research  group. 
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C.  Dynamical  tests 

A  considerable  problem  arises  in  working  with  large,  non-linear  numerical 
models  concerning  the  viability  of  the  final  computer  codes.  That  is,  how  can 
one  feel  confident  that  the  code  is  correctly  performing  the  numerical  integra¬ 
tions  originally  envisioned?  Even  changing  a  working  program  from  machine  to 
machine  frequently  introduces  computational  errors  which  cannot  always  be  de¬ 
tected  by  simple  model  runs.  It  is  necessary,  therefore,  to  subject  such  model 
codes  to  rigorous  testing  procedures  whenever  the  codes  are  modified  or  trans¬ 
ported  to  other  machines.  Such  a  procedure  was  undertaken  and  completed  for  the 
dynamical  portion  of  the  MTM  subsequent  to  introduction  of  the  model  changes 
outlined  in  sections  A  and  B  above.  Similar  checks  were  underway  for  the  ver¬ 
sion  transferred  to  the  AFGL  CDC  6600  at  the  time  of  the  stoppage  of  work  on  the 
model . 

Of  particular  concern  is  the  performance  of  the  non-linear  terms  in  the 
dynamical  sections  of  the  MTM.  We  thus  make  use  of  known  conservative  proper¬ 
ties  of  the  model  to  test  for  ''correctness"  of  solutions  under  various  model 
circumstances.  Some  care,  however,  has  to  be  taken  in  this  procedure  since  it 
is  frequently  very  difficult  to  distinguish  true  model  or  programming  errors 
from  normal  numerical  or  machine  induced  inaccuracies. 

One  series  of  tests  which  have  been  completed  for  the  MTM  involves  running 
the  model  with  the  heating,  frictional  dissipation  and  lower  boundary  vertical 
motion  terms  all  set  to  zero.  Thus,  the  quasi-geostrophic  set  of  dynamical  * 

equations  reduce  to  the  form 
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TJ*  -  -  2n  f  -  Jl*.''2*)  - 

ff  -  -  j(*.n  -  £  v2*  (o.i) 

RV2T  =  V*fV(|J) 

and  we  can  show,  for  example,  that  total  energy  (kinetic  plus  available  poten¬ 
tial)  must  be  preserved  (for  details  of  the  model,  see  the  following  sections). 
Table  0.1  contains  the  results  of  several  runs  under  varying  conditions.  As  a 
base  case.  Run  "A"  was  computed  using  the  normal  N-cycle  scheme  of  Lorenz  (1971) 
with  N  =  4  and  an  internal  time  step  6t  =  1  hour.  We  see  from  the  table  that 
%  0.06%  of  the  initial  model  energy  has  been  lost  after  one  day  of  computations 
and  n,  0.19%  at  the  end  of  two  days.  Thus,  the  energy  has  not  been  preserved 
(which,  of  course,  is  not  unexpected)  and  we  must  ascertain  whether  the  inac¬ 
curacy  is  due  to  our  numerical  approximations  or  results  from  some  more  impor¬ 
tant  physical  or  computational  problem. 

Run  "B"  is  similar  to  "A"  but  we  have  removed  the  non-linear  Jacobian  terms 
from  (0.1).  For  this  case,  the  table  shows  that  the  energy  conserves  much  better 
during  the  first  two  days,  losing  only  n,  0.014%.  From  these  results  it  appears 
that  the  Jacobian  terms  generate  the  major  inaccuracies  in  the  model  runs  but 
it  is  still  not  certain  whether  this  can  be  attributed  to  model  errors  or  to 
numerical  approximations.  One  possibility  would  be  to  change  the  N-cycle  rou¬ 
tine  from  four  to  eight  cycles  per  step  as  an  attempt  to  generate  a  more  accur¬ 
ate  solution.  This  can  be  of  help,  particularly  for  the  linear  parts  of  the 
Jacobian  calculations.  Still  maintaining  6t  =  1  hour  for  the  internal  time 
intervals.  Run  "C"  repeats  the  calculation  of  "A"  for  N  =  8  with  no  improvement 
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Table  0. 1 :  Total  energy  as  a  percent  of  the  initial  total  energy  for  days  0,  1, 
and  2  of  test  Runs  "A"  through  "F".  The  conditions  for  each  run  are  described 
below  the  table. 


Day 

Run  "A"* 

Run  "B"* 

Run  "C"* 

Run  MD"* 

Run  "E"* 

Run  "F"* 

0 

100.000 

100.000 

100.000 

100.000 

100.000 

100.000 

1 

99.943 

99.993 

99.942 

99.999 

99.478 

99.299 

2 

99.812 

99.986 

99.797 

99.996 

99.962 

99.719 

*  All  the  Runs  make  use  of  the  Lorentz  N-cycle  time  stepping  scheme  and,  unless 
otherwise  indicated  below,  the  friction,  heating,  and  lower  boundary  vertical 
motion  terms  are  all  zero.  The  specific  conditions  for  each  run  are: 


Run 

"A": 

Run 

"B": 

Run 

"C": 

Run 

"D" : 

Run 

II  £  M  . 

Run 

II  pil  . 

Uses  the  4-cycle  scheme  with  internal  time  steps  St  =  1  hour. 

Same  as  "A"  but  the  non-linear  Jacobian  (advection)  terms  in  (0.1). 
are  zero. 

Same  as  "B"  but  uses  8-cycles. 

Same  as  "A"  but  St  =  0.2  hours. 

Same  as  "A"  but  the  lower  boundary  vertical  motion  (WB0t)  is 
forced  using  the  results  of  a  Stratospheric  Circulation  Model 
(SCM)  computation. 

Same  as  "E"  but  heating  from  the  SCM  computation  has  been  added. 
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in  the  accuracy  of  the  solutions  (as  seen  in  the  table).  On  the  other  hand, 
when  we  repeat  the  calculation  of  Run  "A"  (4  cycle)  but  with  internal  time  step 
intervals  reduced  to  12  minutes  (6t  =  0.2  hours),  the  accuracy  greatly  improves 
(Run  "0")  with  an  energy  loss  of  only  n,  0.004%  during  the  first  two  days. 
Clearly,  the  small  energy  losses  observed  over  the  first  two  days  of  the  model 
test  runs  are  due  to  numerical  inaccuracies  in  the  time  stepping  scheme  rather 
than  to  coding  errors  in  the  Jacobian  terms. 

To  get  an  idea  of  the  relative  importance  of  the  numerical  errors  detected 
above,  we  ran  two  more  experimental  tests.  The  first  of  these  was  computed 
under  the  conditions  of  Run  "A"  but  with  the  vertical  motion  at  the  lower  bound¬ 
ary  of  the  model  introduced  from  the  results  of  previous  runs  of  the  SCM.  For 
the  second  test  we  added  the  computed  heating  values  from  the  SCM  to  the  lower 
levels  of  the  MTM.  The  results,  shown  as  Runs  "E"  and  "F"  in  Table  0.1,  show 
that  the  energy  changes  introduced  by  these  physical  terms  in  the  model  are  at 
least  as  large  as  the  uncertainties  created  by  the  numerical  procedures  used. 
Thus,  reductions  in  the  time  step  increments  used  for  the  model  to  improve  the 
accuracy  of  the  non-linear  terms  are  not  justified  since  they  would  be  masked 
by  the  forcing  and  boundary  terms. 
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1 .  Basic  dynamical  equations  and  coordinate  system. 

The  horizontal  coordinate  system  will  be  longitude  (positive  eastward)  and 
latitude,  denoted  by  X  and  <f>.  This  dependence  will  be  represented  in  spherical 
surface  harmonics,  except  that  certain  terms,  such  as  part  of  the  heating  and 
photochemistry  will  be  evaluated  point-wise  at  selected  values  of  X  and  <p.  In 
the  vertical  direction  pressure  (p)  will  be  used  as  a  coordinate  with  finite- 
differences  being  employed.  These  pressure  levels  will  be  distributed  at  equal 
intervals  of  log  P  in  order  to  give  roughly  equal  intervals  in  height.  We 
define 

- 

P  =  p  (100  cbar) 

Z  O-l) 

Z  =  -In?,  P  =  e~L  • 

From  the  hydrostatic  relation  dp  =  -pgdz  and  p  =  p/RT,  we  have 

.  <’-2> 

The  vertical  levels  will  be  separated  by  a  uniform  value  of  VZ.  To  the  extent 
that  the  temperature  T  is  approximately  uniform  at  near  surface  values,  a  change 
of  one  in  Z  corresponds  to  a  height  change  of  the  order  of  7  km.  The  bottom  of 
the  atmosphere,  but  not  necessarily  of  the  model,  will,  for  simplicity,  be  taken 
at  Z  =  0,  i.e.,  at  p  =  100  cb  instead  of  at  the  conventional  sea-level  pressure 
of  101.325  cb. 

»  The  dynamical  system  not  only  assumes  hydrostatic  balance,  but  also  a 

"quasi -geostrophic  blance"  in  the  horizontal  equations  of  motion.  Because  we 

% 

must  consider  global  processes  over  the  entire  sphere,  this  balance  must  allow 
for  complete  variability  of  the  Coriolis  parameter  f: 
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(1.3) 


f  =  2ft  s  i  n4> 

ft  =  7.292  x  IQ'5  rad  sec"^ 


The  quasi-geostrophic  balance  in  question  is  obtained  as  follows  (Lorenz, 
Tellus,  1960,  P.  364).  First,  we  divide  the  horizontal  velocity  v  into  a  non- 
divergent  part  k  x  Vip  given  by  a  stream  function  ip  and  a  divergent  part  -Vx, 
given  by  a  velocity  potential  x* 


v  =  k  x  Vip  -  Vx 


(1.4) 


If  the  eastward  and  northward  components  of  v  are  represented  by  u  and  v  and  a 
is  the  radius  of  the  earth,  this  is  equivalent  to 


u  =  a  cos6t—  =  -  1  14 _ J—  lX 

a  cos<^  a  a  cos$  gx 


v  =  a  =  _ I-  14  .  1  1* 

dx  a  cos$  9A  a  d 4> 


(1.5) 


The  vertical  component  of  relative  vorticity,  £,  and  the  horizontal  divergence 


of  v  are  related  to  \p  and  x  by 


t,  =  k  •  curl  v  =  V2\p;  div  v  =  -  v2\ 


where  V2  is  the  horizontal  Laplacian  operator  on  the  sphere. 
The  condition  of  the  quasi-geostrophic  balance  is 


V  •  fVij;  =  gV2z 


(1.6) 


0.7) 


where  g  is  gravity  and  z  is  the  height  of  a  constant  pressure  surface.  [Unless 
noted  otherwise,  all  partial  derivatives  with  respect  to  X,  and  t  (time)  are 
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carried  out  at  constant  pressure  (or  Z)].  The  hydrostatic  relation. 


3z  _  1  _  RT 

93P  “  "  P  "  P 


or 


RT 


(1.8a) 


(1.8b) 


enables  (1.7)  to  be  rewritten  as 

V  •  fv||  =  V2RT  •  (1.9) 

Associated  with  this  relation  (which  is  a  simplified  form  of  the  equation 
obtained  by  taking  the  horizontal  divergence  of  the  equations  of  motion)  is  the 
"vorticity  equation": 

V2||  =  -  k  X  Vf  -  V  (f+V2^)  +  y  •  fVx  +  v  •  (?rxk)  (1.10) 

where  Fr  is  the  horizontal  frictional  force  per  unit  mass. 

The  continuity  equation  (conservation  of  mass)  is 
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The  upper  boundary  condition  at  Z  =  ZtQp  will  be  that  dp/dt  vanishes  there.  Let 
us  define 


X  = 


X  = 


9X 

■  w 


(1.12) 


Equation  (1.10)  can  then  be  rewritten  as 
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V2||  =  -  k  x  Vi)i  •  V(f+V2*)  -  V  •  fV  |pj  +  V  •  (frxk)  .  (1.13) 


If  we  use  Z  =  -InP  as  the  vertical  coordinate,  the  appropriate  vertical 
advection  velocity  is 


1 

P  dt 


(1-14) 


The  continuity  equation  (1.11)  in  terms  of  W  is: 

V-Pv  +  3(PW)/3Z  =  0  •  (1.15) 

From  (1.11),  (1.12)  and  (1.14)  we  get  3[PW  -  V2x]/DP  =  0,  or 

PW  =  V2x  (1.16) 

Boundary  conditions  on  W  are  that  W  vanishes  at  ZtQp  and  that  it  is  given 
from  external  sources  at  the  bottom: 


Z  = 


(1.17) 


Z  =  Zbot:  W  =  Wq ( t , X , )  as  given.  (1.17a) 

Since  Z.  .  is  some  distance  above  the  actual  earth's  surface,  we  must 
bot 

also  specify  the  bottom  dynamical  and  thermodynamical  conditions.  For  this  pur¬ 
pose,  we  will  make  use  of  previous  runs  of  the  model  version  which  includes  the 
surface  as  its  bottom  boundary.  The  results  from  such  a  computation  will  be 
used  to  specify  the  bottom  boundary  temperature  field  (in  space  and  time)  for 
the  present  upper  level  model.  Thus,  we  have 

z  =  Zbot:  T  =  T0(t,A,4>)  as  given.  (1.17b) 
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1 


J  . 


The  bottom  stream-function  field  will  then  be  given  through  the  thermal  wind 
equation. 

Friction  will  be  represented  by  a  vertical  Austausch,  ?r  =  ^  3t/9z  = 
-g3x/3p.  Thus  V*?rxk  =  [V  •  (fj^-ixk)].  We  set  x  =  pK  3(kxVi|j)/3z ,  givinc 


Using  the  "scale  height" 


v.(^xxk)  ■ 

'p  '  pm  3p 

o  o 


RT 

H  =  — 2- 

o  g 


(1.18) 


replacing  p  by  p/RT  and  replacing  g/RT  by  1/H0  we  get 


’  -  F  p  TZ1 
0  0 


To  summarize  the  friction  term  we  can  write 


V-?rxk  =  |p(PF) 


z>0: 

0 


(1.19) 


At  Z  =  Ztop,  F  will  vanish  (no  stress). 

The  next  physical  statement  is  the  thermodynamic  law  d  (entropy)  /dt  = 
rate  of  heating  temperature.  For  our  perfect  gas  system  this  would  be 

cp  3t  [^(Tp’K)]  =  T  ;  K  =  =  7  (K20) 

where  q  is  the  rate  of  heating  per  unit  mass  and  T  the  temperature.  In  terms 


i 
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of  T,  this  becomes 


lb  =  -(^-VX)  •  VT  -  w|J  -  kWT  +  9- 


(1.21) 


We  will,  however,  use  a  simplified  formof  this,  obtained  by  ignoring  ViJj*VT  and 
by  replacing  T  in  W3T/3Z  and  kWT  by  T,  where  T  is  the  horizontal  average: 

T  =  T(p,t)  +  r(X,*,p,t)  / 

i  p/2  r277  _  r  (1-22) 

T  =  cos  d>d<|>  TdX;  T'  =  0 

-tt/2  0  _J 

[This  definition  of  (  )  and  (  ) "  will  be  applicable  to any  variable.]  Th i s 

greatly  simplifies  the  computations,  and  is  reasonably  accurate  because  »  Vx 
and  DT'/DZ  +  kT'  is  generally  small  compared  to  3T/3Z  +  kT.  The  result  is 


3J  . 

3t  " 


kxViJ;  •  VT  -  W(gJ  +  kT)  +  q/cp 


(1.23) 


However,  this  simplification  has  the  result  that  we  can  no  longer  interpret 
(1.23)  as  forecasting  T,  the  horizontally  averaged  T;  this  is  because  the  hori¬ 
zontal  average  of  (1.23)  gives  simply 


whereas  the  horizontal  average  of  the  exact  equation  (1.21)  gives 


3J 

3t 


9-  -  kW'T'  -  1 1-  (P -FT), 


(1.24) 
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showing  the  effect  of  vertical  transports  of  entropy  by  the  motion.  We  expect 
little  change  in  T  from  the  observed  annual  average  T(Z),  however,  either  with 
season  or  with  changes  in  the  ozone  chemistry.  [The  effect  of  the  latter  will 
be  discussed  separately.] 

In  passing,  we  note  that 


<T 


=  RI  fa I  +  g_l 
g  i3z  y 


=  t|-2-  ten(Tp'K)] 


ftj 
R  ' 


(1  25) 

/ 


where  N  is  the  buoyancy  frequency. 

Finally,  we  describe  the  basic  form  of  the  equation  for  the  (number  density) 
mixing  ratio  of  a  trace  substance  such  as  0^.  Define 


where  n^  is  the  number  density  of  the  i-th  trace  substance,  np)  is  the  total  num¬ 
ber  density.  For  levels  below  %  110  km  we  use 


n  =  p/kT 
m 


k  =  Boltzman  constant  = 


1.330x10  ^  kilojoules  deg  ^ 


(1.27) 


Above  'v  110  km,  n(n  =  £n . . 

The  equation  for  dx^/dt  (the  rate  of  change  following  the  motion)  is 
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F 


dt 


3Xi  -  3X; 

jnr +  (kxV,^vx)  •  vxi  +  w^-- 


m 


"dn.j 

,  1  3 

V 

c  P  32 

pKd  3z 

where  (dn^/dt)  is  the  net  rate  of  local  photochemical  generation  of  the  sub¬ 
stance  (number  per  unit  volume  per  unit  time)  and  is  the  vertical  eddy-dif- 

2 

fusion  coefficient  [with  dimensions  length)  f  time].  Kd  will  vary  only  with  P. 

The  vertical  diffusion  term  can  be  rewritten  by  using  the  hydrostatic  equa¬ 
tion  as 


—  TK  Jl  r  3.  r 

3PL  d[RTJ  3P  3Pl' 


K 


4  P^i] 

HJ  3Z  1 


~! 


(1.28) 


where  we  have  again  absorbed  the  variation  of  density  with  T  into  HQ  on  the 
recognition  that  itself  is  not  a  precisely  known  quantity.  Kd  (and  the 
momentum  Austausch  Kj  will  be  prescribed  functions  of  P.  The  equation  for  x 
is  now 


[having  made  use  of  (1.4)  and  (1.15)  to  obtain  the  last  form]. 

The  rate  of  chahge  of  x.j  (the  horizontal  average)  is  obtained  from  the 
horizontal  average  of  (1.30): 


3'X.J 

3t 


dni 


dt 


DP 


[- 


IX1] 

DZJ 


(1.31) 


The  rate  of  change  of  will,  however,  be  obtained  from  a  simplified  form 
of  (1.29),  much  as  was  done  in  the  thermodynamic  equation  (1.23): 


-  3Xi 

at1  =  ■  kxV*  •  -  W3Zl 


i  dni 

r _ 1] 

L"  dt  Jr 


m 


DP 


[- 


n  .32) 


In  contrast  to  f,  where  we  are  for  the  most  part  content  to  take  T  as  given,  we 
must  predict  x\  as  well  as  x(-  Equation  (1-31)  will  therefore  be  used  as  well 
as  (1.32). 

Presumably  (1.33)  need  not  be  applied  every  time  step  in  the  numerical 
integration,  ■<.  being  a  slowly  changing  function  of  time.  However,  the  term 
W'xj  ,nust  be  Put  equal  to  zero  at  P  =  1  to  ensure  no  net  creation  of  xi  by  the 
large  scale  motion. 

A  special  treatment  of  the  minor  species  equation  will  be  necessary  at  cer¬ 
tain  levels.  As  an  example,  lindzen  and  Goody  ( J .  Atmos .  Sc i . ,  1965,  P.  341) 
show  that  the  photodissociation  of  ozone  is  extremely  rapid  at  heights  above  -  45 
km,  with  a  time  constant  becoming  less  than  1  hour.  (They  presumably  use  typi¬ 
cal  values  of  incident  solar  radiation).  The  conventional  methods  of  "time- 
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stepping)  equations  such  as  (1.32)  require  a  computational  time  step  no  longer 
than  the  character) Stic  physical  times  associated  with  terms  on  the  right  side 
of  (1.32).  Since  the  advective  time  scale  is  of  the  order  of  an  hour  or  so,  we 
must  consider  replacing  (1.31)  and  (1.32)  at  these  levels  by  the  equi 1 ibrium 
condition. 

*i  =  (*i)equil  <=>  dF  =  0  (1'33) 
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2 .  Choice  of  vertical  levels. 

We  obtain  equal  intervals  in  Z  =  -In?  (P  =  pressure  •:  100  cb)  by  defining 


Zj  =  AZ(J-j) 


J  =  1,  2,  J. 


(2.1) 


PJ  =  6 


-AZ(J-j) 


j  -  1  is  at  the  "top"  of  our  model  atmosphere,  and  j  =  J  at  the  bottom,  whence 


7  _  Z1  _  Ztop 
Z  '  J^T  "  J^l 


A  convenient  choice  is  obtained  by  choosing 


eAZ  =  r,  r  =  2.12472 
AZ  =  tnr  *  0.753640 


(2.2) 


so  that 


Z1  =  Ztop  = 
P,  =  r-<J-]> 


(2.3) 


Successive  pressure  levels  are  separated  by  (roughly)  6  km  below  the  turbopause. 
The  relations 


pj  =  Pj+1  ■  rPj  (2.4) 

are  useful.  At  these  levels,  the  following  basic  variables  will  be  represented 
j  =1,2,  ...,  J:  Tj,  Wj,  (\ • ) j  toqether  with  the  heating  rate,  the  photochemical 
term,  and  the  vertical  turbulent  fluxes  of  momentum.  At  the  intermediate  levels 


i 


2- 


the  streamfunction  tp.  will  be  represented 

J 


i  =  1  5  1 

j  2»  2’  •  •  •  j  w  2 


For  convenience  in  notation,  however,  tp  will  be  labeled  with  an  interger  sub¬ 
script  according  to  the  convention 

^(P  =  Pj+l/2)  ;  • 

This  results  in  the  scheme  as  seen  in  Figure  2.1. 

Table  2.1  lists  the  values  of  the  more  basic  variables  for  the  choice 
r  =  2.12472,  J  =  26.  Values  of  T  were  taken  from  the  U.S.  Standard  Atmosphere, 

1976  (NOAA,  NASA,  and  USAF).  The  static  stability  parameter  S  ts  defined  later 
in  equation  (3.20). 


r 


Figure  2.1:  Vertical  levels  of  the  model  and  the  location  on  these 
levels  of  the  model  variables. 
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TABLE  2.1 


Pressure,  temperature,  approximate  height, 
and  static  stability  for  model  levels. 


Level 

Z 

(=  -£n(p/1 OOOmb) ) 

p(mb) 

T(°k) 

z(km) 

dZ  c  n  ; 

P 

1 

24.92 

0.15  (-7) 

995.5 

398.0 

289.01 

2 

24.17 

0.32  (-7) 

991  .0 

354.3 

292.70 

3 

23.42 

0.68  (-7) 

981.0 

313.6 

299.13 

4 

22.66 

0.14  (-6) 

962.5 

275.2 

308.45 

5 

21 .91 

0.31  (-6) 

930.5 

240.7 

321.53 

6 

21.15 

0.65  (-6) 

878.5 

210.0 

336.53 

7 

20.40 

0.14  (-5) 

801.5 

183.1 

346.05 

8 

19.65 

0.29  (-5) 

702.0 

161  .0 

345.16 

9 

18.89 

0.62  (-5) 

583.5 

143.0 

327.56 

10 

18.14 

0.13  (-4) 

459.5 

129.0 

288.16 

11 

17.39 

0.28  (-4) 

347.0 

118.9 

233.47 

12 

16.63 

0.60  (-4) 

257.0 

111.4 

162.65 

13 

15.88 

0.13  (-3) 

212.5 

106.0 

100.51 

14 

15.13 

0.27  (-3) 

197.0 

101.0 

71  .20 

15 

14.37 

0.57  (-3) 

190.0 

96.6 

60.91 

16 

13.62 

0.12  (-2) 

187.0 

92.3 

55.41 

17 

12.86 

0.26  (-2) 

187.0 

88.0 

50.76 

18 

12.11 

0.55  (-2) 

191 .0 

83.9 

45.93 

19 

11.36 

0.01 

200.0 

79.5 

45.52 

20 

10.60 

0.02 

208.5 

74.8 

46.62 

21 

9.85 

0.05 

219.5 

69.8 

47.77 

22 

9.10 

0.11 

231 .0 

64.7 

49.07 

23 

8.34 

0.24 

245.0 

59.3 

50.08 

24 

7.59 

0.51 

261  .0 

53.7 

59.96 

1C 

U  J 

6.84 

1  .07 

267.0 

47.8 

80.58 

26 

6.08 

2.28 

254.5 

41  .9 

88.62 

Note:  Levels  21-26  (between  the  dashed  lines)  are  levels  which  overlap  the 
Stratospheric  Circulation  Model. 
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3 .  Non-dimensional  finite-difference  equations 

In  this  section  we  write  the  basic  equations  in  a  non-dimensional  form  (pri¬ 
marily  to  simplify  the  dynamical  computations)  and  simultaneously  introduce  the 
vertical  finite-difference  representation  defined  in  Section  2.  We  define 


U  =  s  i  nf+> 

V(dim)  =  j  7(non-dim) 

V2(dim)  =  V2(non-dim) 

a 

vp ( d i m )  =  2 Qa2  ^(non-dim) 

X(dim)  =  2i'„a:  X(non'dim) 
t(dim)  =  ££  t(non-dim) 

W{dim)  =  2f.W(non-dim) 

T(dim)  =  (4c2a2/R)  T  (non-dim)  +  (4Q2a2/R)  T(non-dim) 


(3.1) 


In  the  last  expression  T  (dim)  is  the  "total"  temperature  in  absolute  degree,, 

T  =  T(Z)  is  the  "standard  atmosphere"  temperature  (also  in  degrees)  given  in  the 
table  at  the  end  of  Section  2,  while  the  cuantity  (4ft2a2/R)  T  (non-dim)  is  the 
(deviation  from  the  horizontal  mean)  variable  T  appearing  in  (1.23),  having  a 
zero  horizontal  average.  [The  total  T  (dim)  is,  of  course,  used  in  all  chemical 
computations] . 


12  =  2i/8. 64x10^  rad  sec 
a  =  6.371x10^  meters 
R  =  287  kd  ton  ^  deg'^ 
Cp  -  ( 7 / 2 )  R 


(3.2) 
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One  day,  (2tt /tt)  secs,  corresponds  to 


At(non-dim)  =  2^(^y)  =  4tt 


The  non-dimensional  72  operator  is 

V2(  )  =  — L,_  J <L  [cosa h 

1  ;  cos 2  j>  IX7^  cos?  a$  o<j)  J 


The  relation 


PW  =  72X 


(3.3) 


(3.4) 


r.i6) 


between  W  and  X  can  be  used  to  eliminate  X  in  favor  of  W  [in  equation  (1.13) 
by  defining  the  inverse  Laplaci an  operator 


We  also  have 


L  =  V2 
X  =  P/.W 


(3.5) 


^  =  V2ip  ;  =  Lt; 


(3.6) 


A  further  convenient  arrangement  is  useful  for  evaluating  terms  of  the  form 
3(PF)/DP,  which  appears  in  the  vertical  diffusion  terms  for  vorticity  and  trace 
substances  and  in  the  term 


in  the  vorticity  equation  (1.13).  We  have 
P 


[,-p(PF)]J 


j±!Z 


9  Fi±Jil 


p 

j+l/2 


'  Pj~V2  Fj-V2 


_  p . 

j-1/2 


(FT)Fj+l/2  '  (FT)Fj-l/2  (3-7) 


where  we  have  made  use  of  (2.4). 
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The  horizontal  advection  of  a  quantity  F  can  be  written  as  the  Jacobian 


VF  =  -  kxVtj;  •  7F 

=  J (F .^) 


3F  3j£  _  3£  9F 
3  A  3y  3  A  3y 


(3.8) 


The  non-dimensional  form  of  the  vorticity  equation  (1.13),  with  regard  to 
the  subscript  labelling  defined  in  Section  2,  together  with  equation  (1.19)  and 
(3.5)  -  (3.8)  is  as  follows: 


For  j  =  1 ,  2,  ....  J-l : 

S1  *  J * 


+  (?^r)Fj+i  •  (?TT)Fj 

(3.9) 

*J  * 

Lcj 

(3.10) 

F!  * 

0 

(3.11) 

FJ  ’ 

■^j-i 

(3.12) 

Fj  ‘ 

Ej(r'j"cj-i)  (j  =  2>  3*  •••» J_1) 

(3.13) 

Eo  3 

(Km)j  v  [H*2nAZ] 

(3.14) 

D  = 

kD  f 

(3.15) 

W1  = 

0 

(3.16) 

WJ  = 

(3.17) 

The  non-dimensional  form  of  the  "thermal  wind  equation"  (1.9)  becomes  for 
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j  -  2 ,  3,  J-l : 

V*nV(,,'ij-'t»j_i )  =  -  V ;>  T  j  A  Z  (3.18) 

The  non-dimensional  form  of  the  thermal  equation  (1.23)  becomes  for 


j  -  2,  3,  J-l: 


2  J!V  Vv0-1> 


where 


SA  +  [C.3ll’a 


R  n 
^]qj 


( 


R 

4TF 


fdl 
a-;  LdZ 


(3.19) 


(3.20) 


is  tabulated  at  the  end  of  Section  2.  Note  that  q.,  the  rate  of  heating  per 

J 

unit  mass,  is  still  in  dimensional  form  in  (3.19). 


The  trace  substance  is,  for 


j  =  j0,  jQ+l,  J-l: 


ST  =  \  J(xj.  '  wj(al>  *  <Fm-)Gj  ■  (Frr,Gj-i  + 


'2ft  Ln  'dt;cJj 
m  0 


°j  ’  °j(xj+r<j>  ;  for  J-J0 . J-2 

°j -  (Kd>j*i/.! :  <2::h;sz> 


(3.21) 


(3.22) 


[The  vertical  diffusion  coefficient  Kd  is  defined  at  the  Z^-levels  corresponding 
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to  j  =  integer  plus  1/2,  whereas  the  vertical  exchange  coefficient  K  for 
ticity,  appearing  in  (3.14),  is  defined  at  interger  values  of  j.] 
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We  define  spectral  solutions  at  arbitrary  level  j  in  the  form 


♦j  ■  £VjVa(A'|!) 
=  K.-jYa(X’“> 

a 

w,  =  yw  ,.y  (x,p) 

J  £  a’j  a 

Tj  - 

"j  -  IVjVa(X-“> 


(4.1) 


and  for  the  trace  substance  equation 


Xj  '  (X.u) 

a 

Gj  ■  i'.’iV1-"' 


(4.2) 


In  terms  of  longitude  (X)  and  latitude  (p),  we  have  defined  members  of  the 
complete  set  of  orthogonal  spherical  harmonics  in  (4.1)  and  (4.2)  using 


it  X 

YjA.y)  =  e  a  Pjy) 


a 


a 


(4.3) 


with 


a  =  n  +  it  (4.4) 

a  a  ' 

denoting  a  vector  index  of  planetary  wave  number  t^  and  degree  n^.  The  Pa(i‘) 
are  Legendre  polynomials  of  rank  and  degree  given  by  u. 


Normal  i^.ition  of  the 


spherical  harmonics  is  such  that  integration  over  the  unit  spherical  surface 
(s)  yields  the  orthogonal  property 


Y  Y*  ds  =  4^6  , , 
,  .  a  v>  a  u> 


(4.5) 


Complex  conjugate  values  are  denoted  by  an  asterisk.  Another  useful  property 
of  the  set  of  spherical  harmonics  is  that  they  satisfy  the  differential  equation 


V‘  Y 

a 


-c  Y  ; 

u  a 


c 

a 


n  (n  +1 ) 
a  a 


(4.6) 


The  complete  set  of  orthonormal  Legendre  polynomials  as  used  in  (4.3)  are  de¬ 
fined  such  that 


P  *  =  P 
a  a 


(4.7) 


and  ail  P  have  been  normalized  such  that 
a 


+1 

-1 


P  P. 

a  3 


dp 


(4.8) 


We  now  want  to  substitute  solutions  (4.1)  and  (4.2)  into  the  non-dimensional 

forms  of  our  model  equations,  multiply  through  with  a  member  of  the  orthogonal 
★ 

set  (say,  Y.  ) ,  and  integrate  the  resulting  relationships  over  the  unit  sphere. 
Application  of  this  procedure  to  the  vorticity  equation  (3.9),  for  example, 
yields  the  desirec  spectral  form  of  this  eo nation. 
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dr 


=  i  -  K  i +  -  i+i 

u  i  ( *J  » J  C  L  '  i  ”1'  -»J+  ■  " 


■  (-Mw 

v-i7  Y-e,J 


E 


r 


Y+c 


^r-l'^V--;  ,j  +  l 


(4.9) 


-(CT>V,j|  +  'r^yd*!  -  <rVl>FvJ 


in  which,  over  the  unit  spherical  surface  s. 


dr 

Y  »J 
dt 


it  :p  • 
i  Y.J 


»?  * 


y,j(*j.“,vv*d» 


>) . 
'3 


‘,nj  S'A  y 


Y*ds 


AY,j  =  4"l  J^j*‘'j)YYdS  ^See  APPendix  A) 


-J-  W 


1  f 


CY-c  Y-c.J  cy+c  Y+c,j 


4-n‘j 


:  [7*u7L(W.)]Y*ds  (See 


Y 


Appendix  B) 


Y  >  J  4-rj 


F . Y*ds 
sJ  Y 


(4.10) 


Similarly,  the  thermodynamic  energy  equation  (3.19),  the  trace  substance  equa¬ 
tion  (3.21),  and  the  thermal  wind  relationship  (3.18)  reduce  to  the  spectral 
forms 
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dT 


-r^’J  =  -B  •  -  S.w  •  + 
dt  y,J  J  y,j 


Cpdfl7 a' 


Y  t  J 


dX,. 


=  -B(^  -  (^)W  -  +  (-iV)G  • 
dt  Y.J  dZ'  Y.J  r-r  y.J 


1 


(“fK 


1 


[, 

-(—) 
n  ^dt'c 

2i.i  i 

S  | 

L  m  J 

Y*ds 
:  Y 


AZ  C  T  .  =  -D  ( 0  *  ,  -tl;  . )  +  E  { tl)  .  . -ib  •) 

Y  Y.J  Y^Y-e.J-1  'Y+e,j; 

where,  for  example. 


?Vj.i 


r  3T  i 


dt  4rr 


3 1 


Y*ds 


%  j 

s/0  = 


dXi 


.dt  y 


Y^ds 


CYTY,J  4ttJ 


(-7:T.)Y*ds 
J  Y 


B 


1 


Y  >  J  8tt 

tx)  =  i_f  ,,u 


,Tj)Y*ds  (See  Appendix  A) 


Y.J  =  8ij  ,Xj)Y*dS  ^See  Appendix  A) 

Vy-cJ  ‘  Vv+r.j  =  '  y  t7*^j^YYds  (See  APPendix  B) 


(4.11; 


3 


(4.12 


In  addition,  we  want  to  determine  the  spectral  form  of  (1.6)  relating  the  verti- 
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cal  component  of  relative  vorticity  (;)  and  the  streamfunction  (y).  It  can  he 
shown  that 


S 


Y.J 


~cy'\  » J 


(4-13) 


or 


(4.14) 


provided  that  in  (6.14)  we  stipulate  -y?*0+i 0  (i.e.,  c  j*0). 

The  spectral  relationships  (4.9),  (4.11),  and  (4.13)  [or  (4.14)]  along  with 
definitions  (4.10)  and  (4.12)  form  a  complete  set  of  equations  for  solution. 
However,  it  is  not  convenient  to  attempt  to  integrate  the  model  in  this  form  as 
there  is  no  explicit  relationship  determining  the  vertical  velocity  field  repre¬ 
sented  by  W.  In  order  to  define  W,  we  want  to  alter  the  thermal  wind  relation¬ 
ship  in  (4.11)  This  development  is  contained  in  the  next  section.  Furthermore, 
specification  of  the  truncation  limits  to  be  used  for  series  solutions  (4.1)  and 
(4.2)  have  not  yet  been  established  and  will  be  discussed  in  a  later  section. 


5.  Determination  of  W  in  the  dynamic  equations 

In  order  to  obtain  an  explicit  description  of  the  vertical  motion  fields  in 
our  model  atmosphere,  we  insert  (4.14)  into  the  thermal  wind  equation  of  (4.11) 
and  differentiate  w.r.t.  time  to  get 


dT 


AZ  c 


Y,J  _ 


D 


Y 


y  d  'C 


Y-e 

E 


X 

“y+£ 


’Y-e, j-l 

dl’Y-e,j 

dt 

dt 

dCY+c 

(5.1) 


dt 


dt 


for  all  levels  j  =  2,3,  _ _  J-l.  We  note  that  (5.1)  does  not  apply  for  the 

cases  y  =  0+iO.  Furthermore,  for  notational  purposes,  we  will  stipulate  that  in 
(5.1)  and  all  future  relationships,  terms  which  require  y-c  =  0+  i 0  or  n 
do  not  exist.  This  applies  equally  to  cases  in  which  y+c  is  not  contained  within 
the  specified  model  truncation  limits. 

Let  us  now  define 


3Y.j  "  "  ~  *Y  >  J  ^  *  AY,j-l  +  - 


+  ^r-T^Y.j  "  ^r-l^FY,j+l 


E5  «  +  [? 


r]q. 


'Y.j  ‘  *Y.j  t  LC^8nTaTrjMY , j 


such  that  using  (4.9)  we  can  write 


(5.2) 


5-1 


J-l 


dt 


d.; 

dt 


’Y,J  _ 


=  a 


1 


D 


■  -  7 — ,-T  [W  •  ,  -  (r+1  )W 

Y,J  Tr-})  C^_r.  L  Y-e»J  + 


+  rV:j+i]  +  tftt  r+;  cVcj-i 

-  (r+l)W  .  ] 

Y  L,J  y+ej+l 


(5.3) 


and,  the  thermodynamic  energy  equation  of  (7.11)  reduces  to 


dT 


dt 


Y  *J 


=  b  .  -  S.W  . 
Y  >J  J  Y  >J 


(5.4) 


Inserting  solutions  (5.3)  and  (5.4)  into  (5.1)  has  the  effect  of  eliminating 
the  time  dependence  of  (5.1)  and  at  any  given  time  we  have 


AZ  CYbr,j  ’  AZ  CYSj\j  =  cy_c  Ve.j  '  Cy+e  Ve,j  - 

,  D  D 

-  -Iz^— Y -  [W  ,  •  ,  -  (r+l)W  ,  .  +  rW  ,  ..,]  + 

(r-D  cY_2ccY_e  y-2e,j-l  Y-2s,J  y-2e,j+lJ 


Ir^T) 


fE  D  ED 

Y ~C  Y  +  Y  Y+C 


c  c  c  c  , 

Y-C  Y  Y  Y+c 


[W  •  ,  -  (r+1 )W  .  +  rW  ...]  - 
L  Y.J-1  v  '  Y.J  Y,J  +  1J 


T^TF  c'^Y+2e^WY+2c,j-l  ‘  (r+1 )WY+2e,j  +  rWY+2c,j  +  l-^ 


or,  if  we  define 
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To  prepare  (5.6)  for  inversion  we  want  to  take  note  of  certain  properties 
of  the  equations  in  order  to  >'educo  the  calculation  to  a  finite  set  of  simple 
matrix  solutions.  Inspection  of  (5.6)  shows  that  the  equations  uncouple  accord¬ 
ing  to  planetary  wave  numbers,  l  In  addition,  within  each  planetary  wave  the 

equations  contain  two  independent  sets;  one  of  even  vector  elements  (n  +  l  all 

Y  Y 

even)  and  the  others  of  odd  vector  elements  (n  +  l  all  odd).  Thus,  to  facili- 

Y  Y 

tate  ease  of  notation,  let  us  define  some  new  sets  of  indices  to  be  applied  to 
(5.6)  by  first  denoting  a  maximum  planetary  wave  number,  L,  for  a  given  spectral 
truncation  as 


L 


(5.8) 


so  that  we  can  designate  K  independent  sets  of  matrix  equations  using  index  k 
where 


k  =  1,  3,  ...,  K;  K  =  2  (L+l ) .  (5.9) 

For  a  given  matrix  set  we  will  determine  k  by  designating 

21  +  1  for  even  vector  sets 
Y 

2(1  +  1)  for  odd  vector  sets 
Y 

Furthermore,  within  each  of  the  K  matrix  equation  sets  it  is  useful  to  designate 
an  element  index,  b^,  where 

bk  =  1.  2,  3 . Bk  .  (5.11) 

Thus,  for  a  given  matrix  set  designated  by  the  subscript  k  we  devise  the  bk 
indices  as  follows: 


(5.10) 


f 
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(1)  For  k  odd  (even  vectors)  let 


N,  =  n.  ) 
k  k  max 


(5.12) 


fo 


r  which  we  consider  only  from  the  set  n^  +  l ^  even.  Then  the  value  for 


an 


individual  b^  is  determined  from 


nk  -  h  *  2  . 

bk  l  '  3f.,0 

k 

Nk  *  ek  +  2 

Bk  2  '  Sl  ,0 

k 


(5.13) 


where  we  ignore  values  of  outside  the  range  indicated  in  (5.11);  i.e.,  when 
k  =  1,  n^  =  0,  and  -  0  we  do  not  include  the  value  b-j  =  0  which  designates 
the  nonallowable  equation  of  (5.1)  in  which  y  =  O+iO  [see  comments  following 
(5.1)]. 

Simi lari ly , 

(2)  For  k  even  (odd  vectors)  let 


^k  nk^max 


(5.14) 


in  which  here  we  consider  only  n^  from  the  set  n^  +  odd.  Then,  we  have 


"k  -  *k  +  1 


Bk  " 


Nk  '  *k  +  1 


(5.15) 


At  this  point  we  want  to  note  an  additional  property  inherent  in  the  spec¬ 
tral  VJ-equations  represented  by  (5.<>).  That  is,  from  definitions  contained  in 


5-5 


(5.5)  and  Appendix  B  we  can  show  that  for  any  given  k, 

E  E  4. 
f (3)  s  Yk 

bk  Sk  Sjc  S{2c 


V 

W 


(5.16) 


We  are  now  prepared  to  convert  (5.6)  to  matrix  form.  To  do  this  we  first 
define  tridiagonal  matrices  V k  as 


:(2)  A3) 


:  ( 3 )  A  2)  A3) 


0  f 


*4,(3) 


V1 


*f(2) 

bk 


’  f  ( 3) 
bk 


f(3)  f (2) 

V1  Bk 


(5.17) 


where  vie  have  made  use  of  (5.8)  -  (5.16).  We  note  from  (5.17)  that  not  only  is 
each  Dk  tridiagonal  but  it  is  also  symmetric.  In  addition,  it  can  be  shown  that 
every  principle  minor  determinate  of  V k  is  positive  and  thus  £>k  can  be  said  to 
be  positive  definite.  These  properties  will  be  discussed  in  more  detail  below. 

To  complete  the  conversion  of  (5.6)  to  matrix  form  we  define  vectors 
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j  ■ 

^lj 

M;j 

‘2,j 

k.j 

’  Rk,j  = 

'bk,j 

k 

.'V 

such  that  (5.6)  can  be  written  in  the  matrix  form 


(5.18) 


Vk.j-l  -  ("’IV'k.j 


rV.  W.  . -  o  -VI,  .  =  R. 
k  k ,  j  +  1  j  k,j  k,j 


j  =  2,  3,  4,  J-1  for  each  k  =  1,  2,  3 . K 


(5.19) 


We  wish  tc  modify  (5.19)  through  diagonal ization  of  each  9^.  However, 
since  eacn  tridiagonal  is  real ,  symmetric  and  po_sji_ti  ve _jiefj_ni  to ,  we  know  that 
all  eigenvalues  of  0^  are  real  and  posi five.  Also,  the  sets  of  eigenvectors 
associated  with  these  eigenvalues  are  orthonormal .  Thus,  if  is  an  MxM  matrix, 
there  exists  a  set  of  real  positive  eigenvalues  (A.)  with  p  =  1,  2,  3,  . ..,  M 

K  p 

associated  with  V.  and  M  sets  of  orthonormal  eigenvectors  q  with  s  =  1,  2,  3, 

K  Pi-3 

. ..,  M.  If  we  let  Qk  represent  the  matrix  of  eigenvectors  associated  with  the 
the  set  (A|Pp  and  matrix  l?k,  we  have 


q12  • 

■  qls  • 

qlm 

q21 

q22  * 

*  q2s  * 

q2m 

>1 

V  * 

•  >  * 

'pm 

qml 

* 

q  ,,  . 
m2 

*  qms  • 

.  q 

mm 

(5.20) 
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such  that 


QkQk  =  QkQk  =  1 


(5.21) 


where  I  is  the  unit  matrix  and  (  )  denotes  transposition.  Define 


A 


k 


'<Vi  0 . ? 

?  ''<v2  i 

I  ; 


(5.22) 


where  then  we  know 


and 


pkQk  QkAk 


Qk^k^k  =  QkQk\  =  Ak 


(5.23) 


We  now  want  to  expand  the  vector  W.  .  in  (5.19)  in  the  form 

K  ,  J 


Wk,j  =  QkVk,j  ;  Vk,j  =  QkWk,j 


(5.24) 


where  we  note  that  V.  -is  also  a  vector. 

k,J 


Inserting  solutions  (5.24)  into  (5.19)  and  multiplying  through  with  Qk  gives 

VkQkVk,j-l  '  (r+1)WkVk,j  +  rQkPkQkVk,j+l  " 

'  '’j^kQkVk,j  =  V*k,j 


or,  from  (5.23),  we  can  write 
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(5 


Vk.j-i  •  t(rt,)  •V'td  \,1  *  r-kvk.M  *  Vk.j 

Now,  we  know  that  there  exists  an  inverse 


Whh  0 . 

0  l/(Xk)2 


.0 


^■Vp 


(5, 


0 .  l/(Aj 


km 


such  that 


\\  - 1 


(5. 


Thus,  if  we  multiply  (5.25)  through  with  ,  (5.c5)  reduces  to  the  form 


'k.j-1  -  +0/k'lVkij  =  \Vk,j 


(5. 


and 


3,  . . . ,  K  we  have  j 

-  2, 

5k,j  E  '  [(r+1)1  + 

,fjAk 

j  \  QkRk,2  ' 

7k,l 

s  Vvkj 

kRk,J-l 

-  rV 

(for  3-_ j <J-2 ) 


(5. 


Using  (5.29)  ,  (5.28)  transforms  to  the  set 


>4*"  k;  if 


~Sk,2Vk,2  +  r  Vk,3  "  Rk,2 

(for  j  =  2)  | 

! 

Vk,j-1  +  Sk,jVk,j  +  r  Vk,j+1  :  Rk,j 

(3<j<J-2)  | 

>  (5.30) 

Vk,J-2  "  Sk,J-lVk,J-l  "  Rk,J-l 

(for  j  =  J-l ) 

1 

from  (5.24)  and  the  boundary  conditions 

of  (5.7)  we  see  that  in  (5.29) 

O 

II 

>• 

■ 

(5.31) 

Vk,J  =  QkWk,J 

J 

We  see  that  for  each  k  the  system  (5.30)  is  tridiagonal  in  j  and  thus  submits 
readily  to  solution  provided  certain  provisions  are  met  (see  Appendix  C  for 
details).  Briefly,  to  invert  (5.30)  we  first  define 


=  c-1 

k,2  '  k  ,2 

(for  j  =  2) 

kj  5  ^Sk,j  "  r  uk  ,j-l  ^ 

(for  3< j<J -1 ) 

k,j  H  "  r  uk,j 

(2<j<J-l) 

(5.32) 


(for  j  =  2) 

yk,j  ’  Vj  (l!k,j  •  yk,j-l>  (for 


•  (5.33) 


6 .  The  model  codes 

The  model  consists  of  four  separate  programs,  three  of  which  are  prelimi¬ 
nary  and  need  to  be  executed  only  once.  The  names  of  the  four  programs  are 
ITC0F1 ,  ITC0F2,  MES0S1 ,  and  MES0S2.  A  brief  description  of  each  of  these  follows. 

ITC0F1 ,  ITC0F2 

ITC0F1  and  ITC0F2  are  used  consecutively  to  generate  and  store  on  the  system 
disk  a  set  of  non-linear  interaction  coefficients  for  use  in  computing  the  non¬ 
linear  Jacobians  in  the  model.  The  definition  and  method  used  for  the  computa¬ 
tion  of  the  interaction  coefficients  are  contained  in  Appendix  A. 

To  run  ITC0F1  ,  use  file  RUNIC1  (see  Figure  6.1).  This  routine  requires  disk 
files  INTC0EF1  and  STRATI  as  input  and  creates  a  file  named  I C 1  OUT  as  output. 

I C 1  OUT  is  used  as  input  by  ITC0F2. 

Ihe  program  RUNIC2  (Figure  6.2)  drives  ITC0F2  and  requires  files  INTC0EF2, 
STRATI,  and  I  Cl  OUT  as  input.  ITC0F2  creates  a  file  IC20UT  on  output  which  con¬ 
tains  the  interaction  coefficient  and  instruction  fields  required  by  the  model. 

MES0S1 

MES0S1  is  an  initializing  program  which  creates  and  stores  all  constants, 
truncation  parameters,  transform  parameters,  and  fixed  fields  required  for  the 
particular  model  configuration  to  be  run.  This  program  must  be  run  prior  to  the 
beginning  of  a  particular  model  experiment  and,  in  general,  calculates  every¬ 
thing  that  can  be  done  for  the  model  in  advance  outside  of  the  main  iterative 
loop. 

MES0S1  is  driven  by  the  file  MFSOSBEG  (Figure  6.3)  and  requires  input  files 
STRATI  and  1C20UT.  On  output,  the  lile  MCSOSW  will  contain  all  of  the  input 


/ 
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fields  retired  by  the  model  except  for  the  model’s  initial  conditions.  In 
addition,  a  set  of  data  which  includes  horizontal  mean  temperatures  and  stabil¬ 
ities  is  required  as  input  as  shown  in  Figure  6.3. 

MES0S2 

MES0S2  represents  the  central  loop  of  the  model.  It  is  called  into  action 
by  the  file  MESOSSTAKT  (Figure  6.4)  which  requires  input  files  STRAT2,  MESOS10, 
INITDEC,  BNDFILE,  and  HEATFILE.  The  last  three  of  these  are  data  files  which 
include  the  initial  conditions  file,  the  lower  boundary  conditions  file,  and 
the  temporary  heating  file  respecti vely.  MESOS10  is  as  described  above  under 
MES0S1. 

On  output,  the  file  MES0UT11  contains  the  complete  history  of  the  model 
run.  The  model  can  be  restarted  as  often  as  required  from  this  history  file  and 
the  subsequent  output  is  appended  onto  the  file.  The  first  records  on  MES0UT11 
contain  fixed  model  parameters.  Subsequent  records  are  each  made  up  of  a  time 
step  number,  the  spectral  vorticity  coefficients,  the  spectral  vertical  velocity 
coefficients  and,  if  applicable,  the  spectral  ozone  mixing  ratio  coefficients 
for  the  particular  time  step. 

An  additional  file  is  provided  in  the  program  to  include  the  spectral  coef¬ 
ficient  values  for  all  the  non-linear  terms  in  the  model  at  each  time  step. 

This  is  assigned  to  "TAPE13"  but  is  currently  not  being  retained  by  the  model  as 
a  permanent  file. 

The  first  pages  of  the  listings  for  the  model  programs  ITC0F1 ,  ITC0F2, 
FiESOSl ,  and  MES0S2  are  shown  in  Figures  6.5  -  6.3  for  purposes  of  reference. 
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Figure  6.7:  hirst  page  of  Program  MES0S1 
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Figure  6.8:  First  page  of  Program  MES0S2. 
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Appendix  A.  Spectral  form  of  Jacobian  terms  and  evaluation  (if  the  a ss oci at e d 
nonlinear  interaction  coefficients . 

Consider  on  the  unit  sphere  the  Jacobian  of  arbitrary  horizontal  global 
scalars  A  and  B  where 


5  fx  If 


3A  3B 
Sp  3A 


(A.  1  ) 


and  A  is  longitude  while  p  is  the  sine  of  latitude.  Expanding  A  and  B  in  terms 
of  spherical  harmonics,  we  have  for  solutions 

A  =  ZaJaU.p). 

a 

8  =  l8aY0(A,U), 

a 

a  =  n  +  i£ 
a  a 

in  which  the  special  properties  of  the  orthonormal  spherical  functions  Y^( \ ,yi) 
are  outlined  in  (4.3)  -  (4.8).  Inserting  solutions  (A. 2)  into  (A.l),  transform¬ 
ing  the  result  to  insure  symmetry  with  respect  to  vector  indices  a  and  g,  and 
writing  in  terms  of  a  single  nonredundant  sun  (for  details  of  these  developments, 
see  Baer  and  Platzman,  1961)  we  arrive  at 


J(A,B)  = 
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for  which  we  define  through  use  of  the  Kronecker  delta. 
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(A. 4) 


The  term 


is  necessary  because  the  two  conjugate  interactions  for  the 


case  11^  =  n^  and  |£  |  =  |£j  ,  assumed  in  the  symmetric  reduction  of  J(A,B)  to 
the  form  of  (A. 3),  are  not  unique  and  ore  of  them  must  be  ignored. 

We  now  multiply  (A. 3)  with  any  arbitrary  member  of  the  orthogonal - 
izing  set,  say  Y*/4ir,  and  integrate  over  the  unit  sphere  to  get 
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and  the  interaction  coefficient,  K  „  is  obtained  from 

Y>B,cx 
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Since  we  intend  to  evaluate  K  „  using  the  "transform"  method  with  integration 
by  exact  Gaussian  quadrature  (see,  for  example,  Eliasen  et  al.,  1970),  a  time 
saving  simplification  can  be  obtained  by  noting  that  the  integral  in  (A. 6)  can 
be  nonzero  only  if  the  integrand  possesses  an  even  parity  with  respect  to  the 
equator.  For  this  condition  we  can  reduce  (A. 6)  to 
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In  order  to  evaluate  (A. 7)  numerically  let  us  define 


A-2 


fo>> 


Qa{\) 


l  P  (y) 

a  av 

dP  (u) 
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(A. 8) 


where  g^  can  be  determined  from  the  Legendre  differential  relationships  in  the 
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Hq  (p)  =  Capp  7T~  -  1  P  7T1 

8, a  8  8  dp  a  a  du 


(A. 10) 


which  can  be  expanded  in  the  form 


H_  (p)  =  Th,  c  Pj.(p) 
3, a  '  £  o,8, a  6'p/ 


(A. 11) 


From  (A. 10)  and  (A. 11)  we  see  that  (A. 7)  can  be  replaced  with 
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However,  if  we  represent  HQ  (y)  at  N  discrete  points  y.  where  k  =  1,  2, 
then  an  exact  quadrature  analog  for  (A. 12)  is  obtained  in  the  form  (see  Eliasen 
et  al. ,  1970) 


Vs.c.  ■ 
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provided 


N  =  (k  +  l)/2 

K  -  ^max  +  2  ^nmax  ”  ^max^  +  2 


(A. 13) 


(A. 14) 


(N  and  k  must  be  intergers)  and  the  latitudes  yk  are  located  at  the  Northern 
Hemisphere  zeroes  of  the  Legendre  polynomial  P^(y)  (including  the  equator  if 
K  is  odd).  In  (A. 13)  the  Wk  represent  the  Gaussian  weights  required  to  maintain 
orthi  jonal ization  of  the  discrete  set  of  Legendre  polynomials  used  in  (A. 13) 
such  that 


j,"kP>k>W  -  Sa,6  <A-,5> 


A  discussion  of  the  evaluation  of  these  Gaussian  weights  is  contained  in  Appen- 
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nd i x  B.  Spectral  representation  of  divergence  terms  of  the  general  form 


In  terms  of  spherical  operators  on  the  unit  sphere  in  which  A  is  longitude 
and  p  is  the  sine  of  latitude  we  have 


V  •  p7A  =  Vp  •  VA  +  pV2A 


=  O-u2)  ^  +  yV2A 


in  which  A  is  an  arbitrary  horizontal  global  scalar  expandable  in  the  form 


A  *  Ta  Y  (A,p) 

L  CL  CL 


Properties  of  the  orthonormal  spherical  functions  Ya(A.p)  are  outlined  in  (4.3) 
(4.8).  Insertion  of  solutions  (B.2)  into  (B.l)  yields 
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where  we  have  defined 
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A  special  case  of  (8.7)  occurs  when  we  consider  scalar  B  in  which 


B  =  V2A 


(B.9) 


where  similar  to  (B.2)  we  can  expand  B  in  the  form 

B  *  Zb  »(»,«). 

a 

Then,  from  (4.6 )>  we  know  that 


and,  in  terms  of  coefficients  b  ,  (B.7)  becomes 
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in  which  we  have  defined 
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(B.13) 


provided  that  in  (B.12)  we  ignore  terms  in  which  c  =  0  (i.e.,  n  =  0). 

y  Y 

Further,  for  both  (B.7)  and  (B.13)  we  must  stipulate  that  all  terms  calling  for 
any  a^,  £,  a^+£,  b^  £  or  b^+£  outside  the  range  of  the  particular  spectral  trun¬ 
cation  chosen  must  also  be  ignored. 
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and  let 


B1  =  C1R1 

By  =  CY(RY-ayBY-1);  2<y<T 


Then,  the  solutions  appear  as 


xr  =  Br 
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provided  all  C  in  (C.5)  are  finite.  That  is,  if 
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We  consider  the  set  of  complete  orthogonal  Legendre  polynomials,  P^(y),  in 

which  t  =  0,  ±1 ,  ±2,  ...  and  n  =  0,  1,  2 .  We  define  this  set,  according 

to  (4.8),  to  be  normalized  such  that 

rl 

P^.(p)dp  =  2«n,n.  (O.D 

-1 


where  y  is  the  sine  of  latitude  or  equivalently,  the  cosine  of  colatitude,  <j>. 
Now  in  order  to  expand  an  arbitrary  function  of  latitude,  say  f(y),  in  terms  of 
the  set  of  Legendre  polynomials  we  let 


f(w)  ■ 

In 


(D.2) 


from  which  the  coefficients,  f^,  are  obtained  through  application  of  (D.l)  such 
that 
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However,  to  be  able  to  transform  at  will  between  spectral  and  grid  point  space, 
it  is  necessary  to  represent  f(u)  at  a  number  of  discrete  points,  y^,  in  which 
k  =  1 ,  2,  3,  . . . ,  N  with  N  being  the  total  number  of  points  lying  within  -l<y<l. 
Thus  at  each  latitude  point,  (0.2)  becomes 

f("k>'n#X>  ■  «>•«> 

•on 

This  means  that  in  order  to  determine  coefficients  fjj  we  must  evaluate  the  inte- 
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grals  in  (D.3)  numerically  and  at  the  same  time  maintain  the  orthogonality  prop¬ 
erties  of  the  discrete  polynomials  representation  in  (D.4).  For  this  purpose, 
integrating  by  quadratures,  we  introduce  a  set  of  Gaussian  weight  functions,  w^, 
such  that 
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and  the  numerical  analog  for  (0.3)  becomes 
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The  remainder  of  this  Appendix  is  devoted  to  the  method  of  evaluation  of  the 
Gaussian  weights,  wk. 

Because  we  know  that  any  given  Legendre  polynomial,  p^(p),  can  be  repre¬ 
sented  by  a  finite  series  in  y  of  at  most  degree  n,  we  can  expand 
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Integrating  (D.8)  by  quadratures  using  (D.5), 
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Equating  (D.8)  and  (0.9)  we  have 
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and  thus  for  any  i  such  that  0  <  i  <  n  +  n'  it  must  hold  that 
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We  see  from  (D.ll)  that  if  we  choose  the  number  of  latitude  points,  N,  such 
that  N-l  =  n+n'  then  utilizing  al 1  i  =  0,  1 ,  2,  . . . ,  n+n  we  can  form  a  set  of 
N  equations  containing  N  unknown  quantities,  w.  ,  for  inversion.  However,  in 
terms  of  colatitude,  <j>,  we  can  show  that  any  cosj^  (j  is  an  integer)  can  be 
expanded  in  the  form 
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Then,  inserting  (D.12)  into  ( D .  1 1 ) , 
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where  we  have  made  use  of  (DJI)  to  eliminate  the  second  term  on  each  side  of 
(D.13).  Again,  as  for  (0.11),  we  see  that  if  we  take  N-l  =  n+n%  we  can  invert 
(D.14)  to  obtain  the  Gaussian  weights. 

As  an  example,  consider  N=3  where  we  select  $^=30°,  <£>2=90°,  and  ^lSO3. 
Then,  from  (D.14)  we  can  construct  the  set  (using  i  =  0,  1,  2) 
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with  solutions 


w1  =  w3  =  4/9 
w2  =  10/9 


(D.  If 


We  note  that  the  solutions  ( D . 16)  are  symmetric  in  about  the  equator.  If  we 
assume  such  symmetry  a  priori  then  ail  equations  in  ( D .  1 4 )  involving  odd  values 
of  i  become  redundant  and  we  can  write  (D.14)  over  the  integration  interval  from 
4>  =  0  to  =  'f/2  as 
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Again,  using  the  example  used  above  in  which  N=3,  ^  =  30°,  and  ^  =  90°,  we 
have  =  2  and  pp  =  1  giving  the  set 


with  solutions 


Furthermore,  if  we  want  to  obtain  w^'s  for  the  entire  pole  to  pole  integration, 
we  need  only  make  use  of  the  symmetry  property 
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which  gives  for  our  example 


(D.20) 


Solutions  (D.20)  are  identical  with  those  of  (D.16). 


